##########################################################################################

library(dplyr)
library(data.table)
library(optparse)
library(MADAM)

##########################################################################################

option_list <- list(
    make_option(c("--report_gene_file"), type = "character") ,
    make_option(c("--score_file"), type = "character") ,
    make_option(c("--bed_file"), type = "character") ,
    make_option(c("--cnv_type"), type = "character") ,
    make_option(c("--type"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20240126_gastric_AFP/analysis/"
    report_gene_file <- "~/20240126_gastric_AFP/analysis/public_ref/SCNA_driver/SCNA_driver_STAD_GENCODE.V19.bed"
    score_file <- "~/20240126_gastric_AFP/analysis/images/SCNA_unpaired/AFP_neg_gene_score.bed"
    bed_file <- "~/20240126_gastric_AFP/analysis/images/SCNA_unpaired/AFP_neg_seg.bed"
    cnv_type <- "Amplification"
    afp_type <- "neg"
    out_path <- "~/20240126_gastric_AFP/analysis/images/SCNA_unpaired"
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

report_gene_file <- opt$report_gene_file
score_file <- opt$score_file
bed_file <- opt$bed_file
cnv_type <- opt$cnv_type
type <- opt$type
out_path <- opt$out_path

dir.create(out_path , recursive = T)
setwd(out_path)

###########################################################################################

gene <- read.delim(report_gene_file,h=F)
met <- read.delim(score_file,h=F)
seg <- read.table(bed_file)

###########################################################################################

## > gene[1:3,]
#     V1        V2        V3 V4              V5        V6              V7  V8
#1  chr1  10686864  10690815  + ENSG00000142655     PEX14  protein_coding AMP
#2  chr1  10696661  10856707  - ENSG00000130940     CASZ1  protein_coding AMP
#3  chr1  11006528  11042094  - ENSG00000175262  C1orf127  protein_coding AMP

#> met[1:3,]
#    V1       V2       V3 V4              V5        V6              V7  V8   V9
#1 chr1 10686864 10690815  + ENSG00000142655     PEX14  protein_coding AMP chr1
#3 chr1 10696661 10856707  - ENSG00000130940     CASZ1  protein_coding AMP chr1
#5 chr1 11006528 11042094  - ENSG00000175262  C1orf127  protein_coding AMP chr1
#       V10      V11      V12      V13 V14
#1 10637416 11197067 0.748591 0.164145 Amp
#3 10637416 11197067 0.748591 0.164145 Amp
#5 10637416 11197067 0.748591 0.164145 Amp

#> seg[1:3,]
#    V1      V2      V3       V4       V5  V6
#1 chr1 3382330 3475479 0.774240 0.166452 Amp
#2 chr1 3475480 3548904 0.782606 0.169442 Amp
#3 chr1 3548905 4169380 0.752494 0.165053 Amp

## met的V12和V13以及seg的V4和V5是log10.q.value.,G.score

## Fisher 的合并方法（Fisher’s method for combining p-values）是一种统计学方法，用于合并多个独立的统计检验的结果。
## 这种方法的主要目的是为了在多重假设检验的背景下，提供一个综合的统计量，以评估这些检验的整体显著性。

###########################################################################################
## amplification ------------------------#
if(cnv_type=="Amplification"){
    gene <- subset(gene, V8 == 'AMP')
    met <- subset(met, V8 == 'AMP' & V14 == 'Amp')
    seg <- subset(seg, V6 == 'Amp')
}else if(cnv_type=="Deletion"){
    gene <- subset(gene, V8 == 'DEL')
    met <- subset(met, V8 == 'DEL' & V14 == 'Del')
    seg <- subset(seg, V6 == 'Del')
}

###########################################################################################
## um-mapped -#
# For genes that did not overlap any GISTIC2.0 segments, the mean G-score of the two neighboring segments were used

umap <- subset(gene, ! V5 %in% met$V5)
umap$log10Q <- NA
umap$G.score <- NA
umap$Type <- NA
umap$chr.amp <- NA
umap$start.amp <- NA
umap$end.amp <- NA
    
for(i in 1:nrow(umap)){
    segs_u <- subset(seg, V1==umap$V1[i] & V2 > umap$V3[i])
    segs_u <- segs_u[which(segs_u$V2 == min(segs_u$V2)),]
    
    segs_l <- subset(seg, V1==umap$V1[i] & V3 < umap$V2[i])
    segs_l <- segs_l[which(segs_l$V3 == max(segs_l$V3)),]
    
    score <- mean(c(segs_l$V5, segs_u$V5))
    p <- fisher.method(matrix(10^(-1*c(segs_l$V4, segs_u$V4)), nrow=1))$p.value
    
    umap$log10Q[i] <- -1*log10(p)
    umap$G.score[i] <- score
    umap$Type[i] <- 'Amp'
    umap$start.amp[i] <- mean(segs_l$V2)
    umap$end.amp[i] <- mean(segs_u$V3)
    
    rm(segs_u, segs_l, score, p, i)
}

umap <- umap[,c('V1','V5','V6','log10Q','G.score','Type','start.amp','end.amp')]
colnames(umap)[1:3] <- c('chr.amp','ID','symbol')

###########################################################################################
## 重复基因
# 取均值
dup <- subset(met,duplicated(met$V5))
dup <- subset(met, V5 %in% dup$V5) 

re <- data.frame()
for(i in unique(dup$V5)){
    sdu <- subset(dup, V5 == i)
    score <- mean(sdu$V13)
    p <- fisher.method(matrix(10^(-1*sdu$V12), nrow=1))$p.value
    
    sre <- data.frame(ID=unique(sdu$V5),
                      symbol = unique(sdu$V6),
                      log10Q = -1*log10(p),
                      G.score = score,
                      Type = 'Amp',
                      chr.amp = unique(sdu$V9),
                      start.amp = min(sdu$V10),
                      end.amp = max(sdu$V11)
                      )
                      
    re <- rbind(re,sre)
    rm(sdu, score, p, i, sre)
}

###########################################################################################
## 非重复
ind <- subset(met, !V5 %in% dup$V5)
ind <- ind[,c('V5','V6','V12','V13','V14','V9','V10','V11')]
colnames(ind) <- c('ID','symbol','log10Q','G.score','Type','chr.amp','start.amp','end.amp')

col_order <- c('chr.amp','ID','symbol','log10Q','G.score','Type','start.amp','end.amp')
umap <- umap[,col_order]
ind <- ind[,col_order]
re <- re[,col_order]

met_am <- rbind(umap, re, ind)
met_am <- met_am[(!is.na(met_am$end.amp)) & (!is.na(met_am$start.amp)) , ]

if(cnv_type=="Amplification"){
    colnames(met_am) <- c('chr.amp','ID','symbol','log10Q','G.score','Type','start.amp','end.amp')
    met_am$Type <- "Amp"
}else if(cnv_type=="Deletion"){
    colnames(met_am) <- c('chr.del','ID','symbol','log10Q','G.score','Type','start.del','end.del')
    met_am$Type <- "Del"
}

write.csv(met_am, file = paste0( "CIN_" , type , "_" , cnv_type , '.Driver_Gscore.csv') , row.names=F,quote=F)